%% Aim
% Perform statiatical analyses based on trial-level pupil data
% dec_pupil: decision pupil

%% Setup
tic;
clc;
clear;
root_path = pwd; 
addpath(genpath(root_path));

nrole=2;
nblock=2;
ntrial=100;
nprice=10;
nlottery=10;

nchoice=2;
choice_label={'Buy','Not Buy'; 'Sell','Not Sell'};

role_label={'Buyer','Seller','Contrast'};
role_color={[226/255,103/255,48/255],[132/255,109/255,177/255],[0.75 0.75 0.75]}; % [Buyer (orange), Seller (purple)]
role_color_light={[236/255,159/255,124/255],[179/255,156/255,192/255],[0.75 0.75 0.75]}; % [Buyer (orange), Seller (purple)]

attribute_label={'Price','Lottery'};
attribute_color={[91/255,155/255,213/255],[112/255,173/255,71/255],[166/255,166/255,166/255]}; % Price vs. Lottery

bias_label={'Val','Res','Val vs. Res'};
bias_color_median={[195/255 205/255 147/255],[192/255 173/255 203/255],[214/255 171/255 128/255],[0.75 0.75 0.75]}; 

%% Load data
sub_data=readtable([root_path,'\data\subject\subject.csv']);
sub_beh_data=readtable([root_path,'\data\subject\subject_behavior.csv']);
sub_ddm_data=readtable([root_path,'\data\subject\subject_ddm.csv']);
beh_data=readtable([root_path,'\data\behavior\behavior.csv']);
pupil_trial_data=readtable([root_path,'\data\pupil\pupil_trial.csv']);
sub_gaze_data=readtable([root_path,'\data\subject\subject_gaze.csv']);
eye_trial_data=readtable([root_path,'\data\eye\eye_trial.csv']);
nsub=height(sub_data);

%% Concatenated data
sub_data = [sub_data, sub_gaze_data(:,4:end)];
trial_data = [beh_data, eye_trial_data(:,9:end), pupil_trial_data(:,9:end)];

%%










%% DecPupil and Choice: Subject-level

for isub=1:nsub
    tdata=trial_data(trial_data.Subject==isub & trial_data.Role==1 & trial_data.ValidPupil==1,:);
    sub_data.DecPupilBuyYes(isub)=nanmean(tdata(tdata.Choice==1,:).DecPupil);
    sub_data.DecPupilBuyNo(isub)=nanmean(tdata(tdata.Choice==0,:).DecPupil);
    
    tdata=trial_data(trial_data.Subject==isub & trial_data.Role==2 & trial_data.ValidPupil==1,:);
    sub_data.DecPupilSellYes(isub)=nanmean(tdata(tdata.Choice==1,:).DecPupil);
    sub_data.DecPupilSellNo(isub)=nanmean(tdata(tdata.Choice==0,:).DecPupil);
end
sub_data.DecPupilBuy=sub_data.DecPupilBuyYes-sub_data.DecPupilBuyNo;
sub_data.DecPupilSell=sub_data.DecPupilSellYes-sub_data.DecPupilSellNo;
sub_data.DecPupilBias=sub_data.DecPupilSell+sub_data.DecPupilBuy;
    
sdata=sub_data;

A=[mean(sdata.DecPupilBuy);median(sdata.DecPupilBuy);std(sdata.DecPupilBuy);min(sdata.DecPupilBuy);max(sdata.DecPupilBuy);std(sdata.DecPupilBuy)/sqrt(height(sdata))];
B=[mean(sdata.DecPupilSell);median(sdata.DecPupilSell);std(sdata.DecPupilSell);min(sdata.DecPupilSell);max(sdata.DecPupilSell);std(sdata.DecPupilSell)/sqrt(height(sdata))];
C=[mean(sdata.DecPupilBuyYes);median(sdata.DecPupilBuyYes);std(sdata.DecPupilBuyYes);min(sdata.DecPupilBuyYes);max(sdata.DecPupilBuyYes);std(sdata.DecPupilBuyYes)/sqrt(height(sdata))];
D=[mean(sdata.DecPupilBuyNo);median(sdata.DecPupilBuyNo);std(sdata.DecPupilBuyNo);min(sdata.DecPupilBuyNo);max(sdata.DecPupilBuyNo);std(sdata.DecPupilBuyNo)/sqrt(height(sdata))];
E=[mean(sdata.DecPupilSellYes);median(sdata.DecPupilSellYes);std(sdata.DecPupilSellYes);min(sdata.DecPupilSellYes);max(sdata.DecPupilSellYes);std(sdata.DecPupilSellYes)/sqrt(height(sdata))];
F=[mean(sdata.DecPupilSellNo);median(sdata.DecPupilSellNo);std(sdata.DecPupilSellNo);min(sdata.DecPupilSellNo);max(sdata.DecPupilSellNo);std(sdata.DecPupilSellNo)/sqrt(height(sdata))];
G=[mean(sdata.DecPupilBuyYes-sdata.DecPupilBuyNo);median(sdata.DecPupilBuyYes-sdata.DecPupilBuyNo);std(sdata.DecPupilBuyYes-sdata.DecPupilBuyNo);min(sdata.DecPupilBuyYes-sdata.DecPupilBuyNo);max(sdata.DecPupilBuyYes-sdata.DecPupilBuyNo);std(sdata.DecPupilBuyYes-sdata.DecPupilBuyNo)/sqrt(height(sdata))];
H=[mean(sdata.DecPupilSellYes-sdata.DecPupilSellNo);median(sdata.DecPupilSellYes-sdata.DecPupilSellNo);std(sdata.DecPupilSellYes-sdata.DecPupilSellNo);min(sdata.DecPupilSellYes-sdata.DecPupilSellNo);max(sdata.DecPupilSellYes-sdata.DecPupilSellNo);std(sdata.DecPupilSellYes-sdata.DecPupilSellNo)/sqrt(height(sdata))];
rowNames = {'mean','median','std','min','max','se'};
colNames = {'Buy','Sell','BuyYes','BuyNo','SellYes','SellNo','BuyDiff','SellDiff'};
dec_pupil.describe = array2table([A,B,C,D,E,F,G,H],'RowNames',rowNames,'VariableNames',colNames);


% one sample t test, one tail
dec_pupil.ottest{1,1}='h';
dec_pupil.ottest{2,1}='p';
dec_pupil.ottest{3,1}='ci';
dec_pupil.ottest{4,1}='stats';
[dec_pupil.ottest{1,2},dec_pupil.ottest{2,2},dec_pupil.ottest{3,2},dec_pupil.ottest{4,2}] = ttest(sdata.DecPupilBuy(:),0);%,'Tail','right');
[dec_pupil.ottest{1,3},dec_pupil.ottest{2,3},dec_pupil.ottest{3,3},dec_pupil.ottest{4,3}] = ttest(sdata.DecPupilSell(:),0);%,'Tail','right');
% paired sample t test, one tail
dec_pupil.pttest{1,1}='h';
dec_pupil.pttest{2,1}='p';
dec_pupil.pttest{3,1}='ci';
dec_pupil.pttest{4,1}='stats';
[dec_pupil.pttest{1,2},dec_pupil.pttest{2,2},dec_pupil.pttest{3,2},dec_pupil.pttest{4,2}] = ttest(sdata.DecPupilBuyYes(:),sdata.DecPupilBuyNo(:));%,'Tail','right');
[dec_pupil.pttest{1,3},dec_pupil.pttest{2,3},dec_pupil.pttest{3,3},dec_pupil.pttest{4,3}] = ttest(sdata.DecPupilSellYes(:),sdata.DecPupilSellNo(:));%,'Tail','right');

disp('Dec pupil');
disp(dec_pupil.describe);
disp('One sample t test, one tail');
disp(dec_pupil.ottest);
disp('t stats buyer');
disp(dec_pupil.ottest{4,2});
disp('t stats seller');
disp(dec_pupil.ottest{4,3});
disp('Paired sample t test, one tail');
disp(dec_pupil.pttest);
disp('t stats buyer');
disp(dec_pupil.pttest{4,2});
disp('t stats seller');
disp(dec_pupil.pttest{4,3});

%% Control Price and Lottery

% data
tdata1=trial_data(trial_data.ValidPupil==1 & trial_data.Role==1,:);
tdata2=trial_data(trial_data.ValidPupil==1 & trial_data.Role==2,:);

tdata1.Price_scaled=rescale(tdata1.Price);
tdata1.Lottery_scaled=rescale(tdata1.Lottery*0.5);
tdata2.Price_scaled=rescale(tdata2.Price);
tdata2.Lottery_scaled=rescale(tdata2.Lottery*0.5);

% Mixed-level regression
% Buyer
glme=fitglme(tdata1,'DecPupil ~ Choice + (Choice | Subject)');
glme=fitglme(tdata1,'DecPupil ~ Price_scaled + (Price_scaled | Subject)');
glme=fitglme(tdata1,'DecPupil ~ Lottery_scaled + (Lottery_scaled | Subject)');
glme=fitglme(tdata1,'DecPupil ~ Choice + Lottery_scaled + Price_scaled + (Lottery_scaled + Price_scaled + Choice | Subject)');
% Seller
glme=fitglme(tdata2,'DecPupil ~ Choice + (Choice | Subject)');
glme=fitglme(tdata2,'DecPupil ~ Price_scaled + (Price_scaled | Subject)');
glme=fitglme(tdata2,'DecPupil ~ Lottery_scaled + (Lottery_scaled | Subject)');
glme=fitglme(tdata2,'DecPupil ~ Choice + Lottery_scaled + Price_scaled + (Lottery_scaled + Price_scaled + Choice | Subject)');

%% Control Price and Lottery, and RT

% data
tdata1=trial_data(trial_data.ValidPupil==1 & trial_data.Role==1,:);
tdata2=trial_data(trial_data.ValidPupil==1 & trial_data.Role==2,:);

tdata1.Price_scaled=rescale(tdata1.Price);
tdata1.Lottery_scaled=rescale(tdata1.Lottery*0.5);
tdata1.RT_scaled=rescale(tdata1.RT);
tdata2.Price_scaled=rescale(tdata2.Price);
tdata2.Lottery_scaled=rescale(tdata2.Lottery*0.5);
tdata2.RT_scaled=rescale(tdata2.RT);

% Mixed-level regression
glme=fitglme(tdata1,'DecPupil ~ RT_scaled + (RT_scaled | Subject)');
glme=fitglme(tdata1,'DecPupil ~ Choice + Lottery_scaled + Price_scaled + RT_scaled + (RT_scaled + Lottery_scaled + Price_scaled + Choice | Subject)');
% extract betas
B=0;
SE=0;
Lower=0;
Upper=0;
for i=1:5
   B(i,1)=glme.Coefficients(i,2);
   SE(i,1)=glme.Coefficients(i,3);
   Lower(i,1)=glme.Coefficients(i,7);
   Upper(i,1)=glme.Coefficients(i,8);   
end
% compare beta: linhyptest
V = glme.CoefficientCovariance;
dfe = glme.DFE;
c = 0;
A = NaN(4,3);
[A(1,1),A(2,1),A(3,1)] = linhyptest(B,V,c,[0 1 -1 0 0],dfe);
[A(1,2),A(2,2),A(3,2)] = linhyptest(B,V,c,[0 1 0 -1 0],dfe);
[A(1,3),A(2,3),A(3,3)] = linhyptest(B,V,c,[0 -1 0 0 1],dfe);
A(4,:)=dfe;
rowNames = {'p','t','r','dfe'};
colNames = {'ChoiceVSPrice','ChoiceVSEVLottery','RTVSChoice'};
disp( array2table(A,'RowNames',rowNames,'VariableNames',colNames));

% Seller
glme=fitglme(tdata2,'DecPupil ~ RT_scaled + (RT_scaled | Subject)');
glme=fitglme(tdata2,'DecPupil ~ Choice + Lottery_scaled + Price_scaled + RT_scaled + (RT_scaled + Lottery_scaled + Price_scaled + Choice | Subject)');
% extract betas
B=0;
SE=0;
Lower=0;
Upper=0;
for i=1:5
   B(i,1)=glme.Coefficients(i,2);
   SE(i,1)=glme.Coefficients(i,3);
   Lower(i,1)=glme.Coefficients(i,7);
   Upper(i,1)=glme.Coefficients(i,8);   
end
% compare beta: linhyptest
V = glme.CoefficientCovariance;
dfe = glme.DFE;
c = 0;
A = NaN(4,3);
[A(1,1),A(2,1),A(3,1)] = linhyptest(B,V,c,[0 1 -1 0 0],dfe);
[A(1,2),A(2,2),A(3,2)] = linhyptest(B,V,c,[0 1 0 -1 0],dfe);
[A(1,3),A(2,3),A(3,3)] = linhyptest(B,V,c,[0 -1 0 0 1],dfe);
A(4,:)=dfe;
rowNames = {'p','t','r','dfe'};
colNames = {'ChoiceVSPrice','ChoiceVSEVLottery','RTVSChoice'};
disp( array2table(A,'RowNames',rowNames,'VariableNames',colNames));


%% DecPupil: Regress on hddm val bias and res bias

sdata=sub_data;

% Normalize data for beta size comparison
sdata.ValBias_scaled=rescale(sdata.ValBias);
sdata.ResBias_scaled=rescale(sdata.ResBias);
sdata.ValBiasLog_scaled=rescale(sdata.ValBiasLog);
sdata.Role1_scaled=rescale(sdata.Role1);

% Main
glm = fitglm(sdata,'DecPupilBias ~ 1 + ValBiasLog_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + ResBias_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + ValBiasLog_scaled + ResBias_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + ValBiasLog_scaled + ResBias_scaled + Role1_scaled');

% Control: val bias as ratio
glm = fitglm(sdata,'DecPupilBias ~ 1 + ValBias_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + ResBias_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + ValBias_scaled + ResBias_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + ValBias_scaled + ResBias_scaled + Role1_scaled');


% Plot
glm = fitglm(sdata,'DecPupilBias ~ 1 + ValBiasLog_scaled + ResBias_scaled ');
disp('Predicted by HDDM parameters');
disp(glm);

B=0;
SE=0;
Lower=0;
Upper=0;
% extract estimates
for i=1:3
   B(i,1)=glm.Coefficients.Estimate(i);
   SE(i,1)=glm.Coefficients.SE(i); 
end
CI=coefCI(glm);
Lower=CI(:,1);
Upper=CI(:,2);
disp('Confidence interval');
disp(CI);

% compare beta: linhyptest
V = glm.CoefficientCovariance;
dfe = glm.DFE;
c = 0;

A = NaN(4,1);
[A(1,1),A(2,1),A(3,1)] = linhyptest(B,V,c,[0 -1 1],dfe);
A(4,:)=dfe;
rowNames = {'p','t','r','dfe'};
colNames = {'ValVSRes'};
temp = array2table(A,'RowNames',rowNames,'VariableNames',colNames);
disp('linhyptest');
disp(temp);

% bar plot
ax=subplot(1,1,1);
set(ax, 'position', [0.15, 0.15, 0.75*1024/1280, 0.8]);
ab=bar([1 2],[B(3,1),B(2,1)],'EdgeColor','none','barwidth',0.5);
xlim([0.2 2.8]);
ab.FaceColor='flat';
ab.CData(1,:) = bias_color_median{1};
ab.CData(2,:) = bias_color_median{2};
hold on
er=errorbar([1 2],[B(3,1),B(2,1)],[SE(3,1),SE(2,1)]);
er.Color = [0 0 0];                            
er.LineStyle = 'none';  
xticklabels({'Val','Res'});
set(gca,'TickDir','out');
xlabel('(scaled)');
ylabel('Beta');
ylim([0 0.6]);
hold off

print(1, '-dtiff', [root_path,'\figure\pupil\dec\ddm\dec_pupil_mean_ddm_regress.tif'], '-r200');
close;

% ratio
B(2)/B(3)

%% Pearson correlation

[r,p]=corr(sdata.DecPupilBias,sdata.vInterceptBias);
[r,p]=corr(sdata.DecPupilBias,sdata.aBias);

%%










%% Plot Gaze bias vs. Pupil bias

sdata=sub_data;

figure('Renderer', 'painters', 'Position', [10 10 370 350]);
ax=subplot(1,1,1);
x=sdata.GazeDurLotteryRatioBiasMinus;
y=sdata.DecPupilBias;
[r,p]=corr(x,y);
scatter(x,y,36,[0.75,0.75,0.75],'filled');
hold on
%ylim([-0.2 0.1]);
%xlim([0 4]);
set(ax,'TickDir','out');
ylabel('Pupillary Bias');
xlabel('Gaze Bias');
box on;

xlim([-0.3 0.3]);
xticks(-0.3:0.1:0.3);
ylim([-0.6 0.8]);
line(xlim,[0 0],'Color',bias_color_median{1},'LineStyle','--');
line([0 0],ylim,'Color',bias_color_median{2},'LineStyle','--');

print(1, '-dtiff', [root_path,'\figure\pupil\gaze_vs_pupil.tif'], '-r200');
hold off;
close;

% save
A=[r; p];
rowNames = {'r';'p'};
colNames = {'Bias'};
disp(array2table(A,'RowNames',rowNames,'VariableNames',colNames));

%% Plot Gaze bias vs. Pupil bias

sdata=sub_data;

figure('Renderer', 'painters', 'Position', [10 10 370 350]);
ax=subplot(1,1,1);
x=sdata.FirstGazeLotteryRatioBiasMinus;
y=sdata.DecPupilSell+sdata.DecPupilBuy;
[r,p]=corr(x,y);
scatter(x,y,36,[0.75,0.75,0.75],'filled');
hold on
%ylim([-0.2 0.1]);
%xlim([0 4]);
set(ax,'TickDir','out');
ylabel('Pupillary Bias');
xlabel('First Gaze Bias');
box on;

xlim([-0.4 0.8]);
ylim([-0.6 0.8]);
line(xlim,[0 0],'Color',bias_color_median{1},'LineStyle','--');
line([0 0],ylim,'Color',bias_color_median{2},'LineStyle','--');

print(1, '-dtiff', [root_path,'\figure\pupil\gaze_first_vs_pupil.tif'], '-r200');
hold off;
close;

% save
A=[r; p];
rowNames = {'r';'p'};
colNames = {'Bias'};
disp(array2table(A,'RowNames',rowNames,'VariableNames',colNames));

%% save
writetable(sub_data,[root_path,'\data\subject\subject_pupil.csv']);

%%
toc;










%% DDM Difficulty

tdata1=trial_data(trial_data.ValidPupil==1 & trial_data.Role==1,:);
tdata2=trial_data(trial_data.ValidPupil==1 & trial_data.Role==2,:);

tdata1.Price_scaled=rescale(tdata1.Price);
tdata1.Lottery_scaled=rescale(tdata1.Lottery*0.5);
tdata1.RT_scaled=rescale(tdata1.RT);
tdata2.Price_scaled=rescale(tdata2.Price);
tdata2.Lottery_scaled=rescale(tdata2.Lottery*0.5);
tdata2.RT_scaled=rescale(tdata2.RT);

tdata1.EVAbs=abs(tdata1.Price-tdata1.Lottery./2);
tdata2.EVAbs=abs(tdata2.Price-tdata2.Lottery./2);

glme=fitglme(tdata1,'DecPupil ~ EVAbs + (EVAbs | Subject)');
glme=fitglme(tdata2,'DecPupil ~ EVAbs + (EVAbs | Subject)');

glme=fitglme(tdata1,'DecPupil ~ Choice + Lottery_scaled + Price_scaled + RT_scaled + EVAbs + (RT_scaled + Lottery_scaled + Price_scaled + Choice + EVAbs | Subject)');
glme=fitglme(tdata2,'DecPupil ~ Choice + Lottery_scaled + Price_scaled + RT_scaled + EVAbs + (RT_scaled + Lottery_scaled + Price_scaled + Choice + EVAbs | Subject)');

%% DDM Difficulty

tdata1=trial_data(trial_data.ValidPupil==1 & trial_data.Role==1,:);
tdata2=trial_data(trial_data.ValidPupil==1 & trial_data.Role==2,:);

tdata1.Price_scaled=rescale(tdata1.Price);
tdata1.Lottery_scaled=rescale(tdata1.Lottery*0.5);
tdata1.RT_scaled=rescale(tdata1.RT);
tdata2.Price_scaled=rescale(tdata2.Price);
tdata2.Lottery_scaled=rescale(tdata2.Lottery*0.5);
tdata2.RT_scaled=rescale(tdata2.RT);

for i=1:nsub
    weight=sub_ddm_data.ValBias(i);
    tdata1.EVAbsDDM=abs(tdata1.Price-weight.*tdata1.Lottery./2);
    tdata2.EVAbsDDM=abs(tdata2.Price-weight.*tdata2.Lottery./2);
end
tdata1.EVAbsDDM_scaled=rescale(tdata1.EVAbsDDM);
tdata2.EVAbsDDM_scaled=rescale(tdata2.EVAbsDDM);

glme=fitglme(tdata1,'DecPupil ~ EVAbsDDM_scaled + (EVAbsDDM_scaled | Subject)');
glme=fitglme(tdata2,'DecPupil ~ EVAbsDDM_scaled + (EVAbsDDM_scaled | Subject)');

glme=fitglme(tdata1,'DecPupil ~ Choice + Lottery_scaled + Price_scaled + RT_scaled + EVAbsDDM_scaled + (RT_scaled + Lottery_scaled + Price_scaled + Choice + EVAbsDDM_scaled | Subject)');
glme=fitglme(tdata2,'DecPupil ~ Choice + Lottery_scaled + Price_scaled + RT_scaled + EVAbsDDM_scaled + (RT_scaled + Lottery_scaled + Price_scaled + Choice + EVAbsDDM_scaled | Subject)');

%%
